Time-dependent complete-active-space self-consistent field method for multielectron dynamics in 

intense laser fields 
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The time-dependent complete-active-space self-consistent-field (TD-CASSCF) method for the description of 
multielectron dynamics in intense laser fields is presented, and a comprehensive description of the method is 
given. It introduces the concept of frozen-core (to model tightly bound electrons with no response to the field), 
dynamical-core (to model electrons tightly bound but responding to the field), and active (fully correlated to de- 
scribe ionizing electrons) orbital subspaces, allowing compact yet accurate representation of ionization dynam- 
ics in many-electron systems. The classification into the subspaces can be done flexibly, according to simulated 
physical situations and desired accuracy, and the multiconfiguration time-dependent Hartree-Fock (MCTDHF) 
approach is included as a special case. To assess its performance, we apply the TD-CASSCF method to the ion- 
ization dynamics of one-dimensional lithium hydride (LiH) and LiH dimer models, and confirm that the present 
method closely reproduces rigorous MCTDHF results if active orbital space is chosen large enough to include 
appreciably ionizing electrons. The TD-CASSCF method will open a way to the first-principle theoretical study 
of intense-field induced ultrafast phenomena in realistic atoms and molecules. 

PACS numbers: 32.80.Rm, 31.15.A-, 42.65.Ky 



I. INTRODUCTION 

The advent of the chirped pulse amplification (CPA) tech- 
nique [ 1 1 has enabled the production of femtosecond laser 
pulses whose focused intensity easily exceeds 10 14 W/cm 
and even reach ~ 10 22 W/cm |2]44). Exposed to visible- 
to-mid-infrared pulses with an intensity typically higher than 
10 14 W/cm , atoms and molecules exhibit nonperturbative 
nonlinear response such as above-threshold ionization (ATI), 
tunneling ionization, high-order harmonic generation (HHG), 
and nonsequential double ionization (NSDI) J3] |6|. HHG, 
for example, represents a highly successful avenue toward 
an ultrashort coherent light source in the extreme-ultraviolet 
(XUV) and soft x-ray regions 0. The development of 
these novel light sources has opened new research possibili- 
ties including ultrafast molecular probing [9-11], attosecond 
science |[T2HT4ll . and XUV nonlinear optics |[T5l[T6l . 

In parallel with the progress in experimental techniques, 
various numerical methods have been developed to explore 
atomic and molecular dynamics in intense laser fields. While 
direct solution of the time-dependent Schrodinger equation 
(TDSE) provides exact description, this method is virtually 
unfeasible for multielectron systems beyond He lfP7H29l and 
H2 (30—32 1 . As a result, single-active electron (SAE) approx- 
imation is widely used, in which only the outermost electron 
is explicitly treated, and the effect of the others, assumed 
to be frozen, is embedded in a model potential. This ap- 
proximation, however, fails to account for multielectron and 
multichannel effects ifTUl [3314371 in high-field phenomena. 
Thus, alternative many-electron methods are required to catch 
up with new experimental possibilities. For example, Cail- 
lat et al. [38 1 have introduced the multiconfiguration time- 
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dependent Hartree-Fock (MCTDHF) approach (see below) to 
study correlated multielectron systems in strong laser fields. 
Although they have presented the results for up to six-electron 
model molecules, its computational time prohibitively in- 
creases with the number of electrons. Another interesting 
route is the time-dependent configuration-interaction singles 
(TD-CIS) method, implemented by Santra and coworkers. 
Il39l l40ll . in which the many-electron wavefunction is ex- 
panded in terms of the Hartree-Fock (HF) ground state and 
singly excited configuration state functions (CSF). The TD- 
CIS method has an advantage to give a clear one-electron pic- 
ture for multichannel ionizations. However, its applications 
are limited to the dynamics dominated by single ionizations, 
with an initial state described correctly by the HF method. 
Time-dependent density functional theory (TDDFT) [41-43), 
though attractive for its low computational cost, delivers only 
the electron density, not the wavefunction, rendering the defi- 
nition of observables difficult. More seriously, it is difficult 
to estimate and systematically improve the accuracy of the 
exchange-correlation potential. 

Among the more recent developments, the orbital-adapted 
time-dependent coupled-cluster (OATDCC) method proposed 
by Kvaal iPHTl is of particular interest, which pioneered the 
time-dependent coupled-cluster (CC) approach with biviria- 
tionally adapted orbital functions and excitation amplitudes. 
The fixed-orbital CC method was also implemented by Huber 
and Klamroth B31 . The time-dependent CC approach should 
be a promising avenue to the time-dependent many-electron 
problems in view of the spectacular success of the stationary 
CC theory. However, it seems to require further theoretical 
sophistications to make a rigorous numerical method on it. 
Hochstuhl and Bonitz proposed the time-dependent restricted- 
active-space configuration-interaction (RASCI) method |46|, 
in which the total wavefunction is expanded with fixed Slater 
determinants compatible to the RAS constraints known in 
quantum chemistry ll47l . The TD-RASCI method, with clev- 
erly devised spatial partitioning, has been successfully applied 



to helium and beryllium atoms [46|. A disadvantage of the 
method is the lack of the size extensivity [48 49 1 discussed 



symbolically as 



in Sec. Ill which is a common problem of truncated CI ap- 



proaches, except the simplest TD-CIS [39|. 

In this work, we propose a flexible ab initio time-dependent 
many-electron method based on the concept of the complete 
active-space self-consistent field (CASSCF) [50-53]. Our ap- 
proach is derived from the first principles, and simultaneously, 
fills the huge gap between the MCTDHF method and SAE ap- 
proximation. 

Most of the ground-state closed-shell systems are described 
qualitatively well by the HF method [4-8 1 . However, the time- 
dependent Hartree-Fock (TDHF) method cannot describe ion- 
ization processes ||54l , since it enforces to keep the initial 
closed-shell structure. Instead, at least two orbital functions 
are required to describe the field ionization of a singlet two- 
electron system. The spatial part of the total wavefunction in 
such approximation reads 

tt(l, 2) <x V>i(l)V> 2 (2) + ^2(1)^1(2) da) 

= Ci&(l)0i(2) + C 2 02(l)&(2), (lb) 

The first form is known as the extended Hartree-Fock (EHF) 
ll55H57l wavefunction, and has been successfully applied to 
the intense-field phenomena for two-electron systems Il55ll57l 
58 1 . The second form is obtained by the canonical orthogonal- 
ization of non-orthogonal orbitals ipi=i,2 [48], and an example 
of the configuration-interaction (CI) wavefunction |48l |49l . It 
is clear that not only the CI coefficients {Ci(i}} but also the 
orbitals {4>i(t)} have to be varied in time in order to prop- 
erly describe the ionization. Thus we need the multiconfig- 
uration Hartree-Fock (MCHF) or the multiconfiguration self- 
consistent field (MCSCF) wavefunction 08] 09) |53l, where 
both the CI coefficients and the shape of the orbitals are the 
variational degrees of freedom. 

This idea has been realized for many-electron systems by 
the MCTDHF method CODED, in which the total wavefunc- 
tion is expanded in terms of Slater determinant (or CSF) bases, 
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where both CI coefficients {Ci(t)} and bases {$/(£)} are 
allowed to vary in time. The Slater determinants are con- 
structed (in the spin-restricted treatment, see Sec. |II A) > from 
2n spin-orbitals {<f> p ;p = 1,2, • • n}<g){a,f3} , where {4> p (t)} 
are time-dependent spatial orbital functions and a (j3) is the 
up- (down-) spin eigenfunction. See also Refs. I160H661 for the 
MCTDHF method, and Beck et al |67 1 and references therein 
for the precedent multiconfiguration time-dependent Hartree 
method for Bosons. 

Despite its naming, which in principle refers to the general 
multiconfiguration wavefunction of the form Eq. |2]) with flex- 
ible choice of expansion bases (range of summation /), pre- 
vious implementations of the MCTDHF method (except the 
fixed-CI formulation of Ref. |66|) were limited to the full-CI 
expansion; the summation / of Eq. |2) is over all the possible 
ways of distributing N electrons among the 2n spin-orbitals. 
More intuitively, we write such an MCTDHF wavefunction 
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which is understood to represent the A-electron full-CI wave- 
function using n time-dependent orbitals. Though powerful, 
the MCTDHF method suffers from severe limitation in the ap- 
plicability to large systems, since the full-CI dimension grows 
factorially with increasing N. 

It is reasonable to expect that in a large molecule interacting 
with high-intensity, long-wavelength lasers, the deeply bound 
electrons remain non-ionized, while only the higher-lying va- 
lence electrons ionize appreciably. For the bound electrons, a 
closed-shell description of the HF type would be acceptable as 
a first approximation |68 1. On the other hand, correlated treat- 
ment is required for ionized electrons to describe the seamless 
transition from the closed-shell-dominant initial state into the 
symmetry-breaking continuum (discussed in Sec. |IIIB] ). 

The CASSCF method [50-53] provides an ideal ansatz 
for such a problem. It introduces the concept of core and 
active orbital subspaces, and spatial orbitals participating in 
Eq. <f3T> are classified into these subspaces. The core orbitals 
are forced to be doubly occupied all the time, while the active 
orbitals are allowed general (0, 1, or 2) occupancies. Thus, 
the CASSCF wavefunction is written symbolically as 
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(4b) 



where factors ( |4a[ > and ( |4b| i represent core- and active- 
subspaces, respectively, and the total wavefunction is properly 
antisymmetrized. See Sec. II A for the rigorous definition. 
The nc core orbitals describe Ac = nc/2 core electrons 
within the closed-shell constraint, while the N\ active elec- 
trons are fully correlated using % active orbitals. Whereas, 
in general, all the orbitals are varied in time, it is also possi- 
ble to further split the core space into frozen-core (fixed) and 
dynamical-core (allowed to vary in time in response to the 
field) subspaces. See Fig. [T] which illustrates the concept of 
the orbital subspacing. 

The equation of motion (EOM) of our method, called TD- 
CASSCF, is derived based on the time-dependent variational 
principle (TDVP) [69-71], as detailed in Sec. [TT] It guaran- 
tees the best possible solution using the total wavefunction ex- 
pressed as Eq. (HJ. The fully correlated active space allows to 
describe ionization processes which involve the strong corre- 
lation due to the breaking closed-shell symmetry (discussed in 
Sec. |IIIB) . It can also include multichannel and multielectron 
effects (discussed in Sec. IIICl. The dynamical core orbitals 
account for the effect of field-induced core polarizations. In 
whole, the TD-CASSCF method enables compact yet accurate 
representation of multielectron dynamics, if the active space 
is chosen correctly according to the physical processes of in- 
terest. 

This paper proceeds as follows. In Sec. [Til the details of 
the TD-CASSCF method are described. Then in Sec.|mJ the 
TD-CASSCF method is numerically assessed for ionization 
dynamics of one-dimensional multielectron models. Finally, 



Sec. IV concludes this work and discusses future prospects. 



Appendix [A] provides the definition and computational de- 
tails for real-space domain-based ionization probabilities. The 
Hartree atomic units are used throughout unless otherwise 
noted. 



II. THEORY 

In this section, we derive the EOMs for the TD-CASSCF 
method. To this end, first we give the rigorous definitions of 
the MCTDHF, MCSCF, CASSCF, and general multiconfig- 
uration ansatz for the total wavefunction in Sec. |HA| Next 
in Sec. |HB| the EOMs of orbitals and CI coefficients for the 
general multiconfiguration wavefunctions are discussed by re- 
viewing the work of Miranda et al ll66l . Then we special- 
ize the general formulation to the TD-CASSCF method in 
Sec. |HC| to derive the explicit EOMs, and discuss its com- 



putational aspects in Sec. II D 



A. Multiconfiguration wavefunctions 

Our formulation is determinant-based l47l l65l within the 
spin-restricted treatment, i.e., using the same spatial orbitals 
for up- and down-spin electrons. We define n occupied spa- 
tial orbitals {4> p ;p = 1, 2, • • -,n} and Ni, — n virtual orbitals 
{<fia, a = n + l,n + 2, • • •, Nt,}, here Nj, is the dimension 
of the spinless one-particle Hilbert space, determined by e.g., 
the number of grid points to discretize the orbitals or the num- 
ber of intrinsic single-particle basis functions to expand the 
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FIG. 1. Pictorial explanation of the TD-CASSCF concept, illustrat- 
ing a 12-electron system with two frozen-core (orbitals 1 and 2), two 
dynamical-core (3 and 4), and four active orbitals (5 to 8). Solv- 
ing the orbital-EOM guarantees the variational splitting, in a time- 
dependent sense, of core-active, core-virtual, and active-virtual or- 
bital subspaces. 



orbitals. The indices {p, q.r, s}, {a, b}, and {/i, is, A, 7, 8} 
are used to label occupied, virtual, and general (occupied + 
virtual) orbitals, respectively, and {er, r} € {a, f3} label spin 
eigenfunctions. In the followings, Einsteins's summation con- 
vention is applied for repeated upper and lower orbital indices 
within a term, with summation ranges implicit as above. The 
orbitals are assumed to be orthonormal all the time, 



(Jt)\M*))= dr<FJt)Mt) = %, 



(5) 



where 8% denotes the Kronecker delta. The Slater determi- 
nants $i(t) in Eq. (pi are constructed, as usual, from the oc- 
cupied orbitals {0 p }7and given in the occupation number rep- 
resentation 149, 651 as 
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where alg. (and a^ a appearing below) is the Fermion creation 
(destruction) operator for 2iV& spin-orbitals {</> M } <8> {a,/?}. 
The integer array {i? = 0,1} specifies the occupancy of 
each spin-orbital, where X)«=i K = N a , with N a being the 
number of er-spin electrons, and N = N a + N 13 . The time- 
dependence of the occupation number vector \I(t)) is implicit 
in the spatial orbitals. 

We focus on the dynamics induced by spin-independent ex- 
ternal fields, and the initial wavefunction is assumed to be the 
spin eigenfunction. Therefore the total and projected spin op- 
erators, S 2 , S z are the constants of motion. Each Slater deter- 
minant is the eigenfunction of S z with eigenvalue N a — N 13 , 
while not generally of S 2 . The total wavefunction is spin- 
adapted, however, since the initial state is prepared by the 
variational optimization of CI coefficients, which automati- 
cally gives proper spin combinations. 

Throughout this paper, the term MCTDHF is used for the 
method based on the full-CI expansion using n occupied or- 
bitals: 



I* 



ripci 

MCTDHF (*)) = J2 C'WW*)}' 
I 



(7) 



with I varying freely in the full-CI space IIfci, spanned by 
all the determinants generated from the 2n occupied spin- 
orbitals. This is the current standard of the MCTDHF method 
as mentioned in Sec. [I] Note that Miranda et al [66| used the 
term MCTDHF in a broader sense to denote approaches based 
on the general multiconfiguration wavefunctions. To be defi- 
nite, we call the general ansatz as MCSCF: 



|*MCSCF(i))=VC / (i)|/(i)) 



li 



(8) 



with the general CI space n defined as any arbitrary subspace 
of n FC i, n = {\I) G II F ci; Cj(t) =£ 0}. A trivial exam- 
ple of this class is the single-determinant HF wavefunction 
for closed-shell singlet or open-shell high-spin states. The 
only nontrivial applications of Eq. dSl to the time-dependent 



problems made thus far is the general open-shell TDHF ap- 
proaches formulated in Ref. |66|, in which the CI coefficients 
are determined by the spin-symmetry and time-independent. 

The most successful MCSCF method in quantum chemistry 
is the CASSCF (also known as fully optimized reaction space) 
method introduced in Sec.|I][Eq. Q], in which the CI expan- 
sion is limited to the space spanned by Slater determinants that 
include nc doubly occupied core orbitals, called the CASCI 
space IIcas: 



I* 



CASSCF 
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Ci(t)\I(t)), 
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a\Ap 



n («L) /r (< 



IS' 



(9) 



(10) 



(=i 



t=n c +l 



with J2 t 1° = N A being the number of cr-spin active elec- 
trons satisfying 7V£ = N a - N c /2 and N A = N% + JVf. 
Hereafter, we use orbital indices {i,j,k,l} for core- and 
{t, u, v, w, x} for active-orbitals, while keeping {p, q, r, s} for 
general occupied (core + active) orbitals. Following the con- 
vention in the electronic structure theory IJ481 |49l , we use 
the acronym CASSCF(7V A , n A ) to denote the CASSCF wave- 
function with Na active electrons and ua active orbitals. The 
MCTDHF wavefunction with n occupied orbitals is identical 
to CASSCF(N, n) and denoted as MCTDHF(n). See Fig. [T] 
Eqs. ((3) and Q in Sec. [I] and Eqs. (|49]> and ^ in Sec.[H^for 
intuitive understanding of these notations. 



E% = J2 «^«L«7ra^ = EgE. 






(15) 



Following the TDVP I69H7T1 . the action integral S, 



1 dt(V\ ( H - i^- ) \^l> 



at 



is made stationary, 

SS= f ' dtUs^\ (h\^) - i\i>)) 

Jt 



(16) 



'MH + i 



= 0, 



(17) 



with respect to allowed variations S"i> of the total wavefunc- 
tion, where 'J = d^/dt. The time derivative of 6*f> is in- 
tegrated out by part, assuming S^f(to) — 6^(ti) — 0. See 
Ref. 11711 for the formal discussion on the validity of this pro- 
cedure. By taking the orbital orthonormality into account, the 
variations and the time derivatives of an orbital <f>^ can be writ- 
ten as 66, , 



b^A" and ic 



f> v Ru> respectively, with 



R£ = i{4>u,\<j>i>}- The matrix A is anti-Hermitian, while R 
is Hermitian [66 1 . Then, the allowed variation and the time 
derivative of the total wavefunction are compactly given by 



|,y*) = 2li)*Cj + A|*>, 



(18) 



B. Equations of motion for MCSCF wavefunctions 

Recently, Miranda et al discussed EOMs for MCSCF wave- 
functions 1 66 1 . Although their main motivation was the fixed- 
CI formulations, they also presented important equations ap- 
plicable to the general MCSCF wavefunctions (See Sec. IV of 
Ref. 1 66]). Here we follow the essentials of their development 



to obtain Eqs. (20 1 and (21 1 below. 

The spin-free second-quantized Hamiltonian is given by 



1 



h = Km + ~ggEfi 



(ii) 



where b% and gfy are the one- and two-electron Hamiltonian 
matrix elements, 



w,= 



^*(r)ft(r,V r )^(r), 



(12) 






dridr 2 0* (r 1 )6* x (r 2 )V ee (r 1 , r 2 ) 

x6 u (ri)6-y(r 2 ), (13) 



with h consisting of kinetic, nucleus-electron, and external 
laser terms, and V ee being the electron-electron interaction, 
and 



^ = E« 



(14) 



i\y) = ij2\i)c I + m, 



(19) 



where 6C1 and Cj are the variation and the time derivative of 
CI coefficient C/, and A = A£J5£, R = R^E^. Inserting 
Eqs. (19 1 and ( fl9]> and their Hermitian conjugates (<J$| and 
— i($| into Eq. (|17[), and requiring the equality for individual 
variations <5Cj (tjand A|(t), after some algebraic manipula- 
tions 1661 we have 



iC I = (I\Sm, 



(20) 



(V\H(1 - il)E^) - (*|££(1 - fl)H\t>) = 0, (21) 

where H = H — R, and II = J2ien 1-0(^1 I s tne configu- 
ration projector onto the gen eral CI space II. The system of 
equations, Eqs. (20 1 and (21 1, is to be solved for iCj and R, 
which determine the time-dependence of CI coefficients and 
orbitals, respectively. In Ref. [66|, these equations appeared 
as an intermediate to derive the MCTDHF equation, rather 



than as the final result. Here we emphasize that Eqs. (20 1 and 



pT} are valid for general MCSCF wavefunction fit into the 
form of Eq. (|8j. Equation pi) is also extensively discussed 
by Miyagi and Madsen in their recent development of MCT- 
DHF method with restricted CI expansions |72|. 



C. TD-CASSCF equations of motion 



1. Orbital equations of motion 



Now we apply the CASSCF constraint defined in Sec. 



to the general orbital-EOM derived in Sec. II B Equation ( 



ITA| 
2TJ7 

with IT replaced by lie as > reduces to a trivial identity for an 
orbital pair {fi, v} belonging to a same orbital subspace (core, 
active, or virtual), since the singly replaced determinants, 
Ejj\l) = 28) 1 1) , El 1 1) , or Eg \ I) = 0, either fall within n c as 

or vanish, and the configuration projector 1— IIcas eliminates 
such contributions. We refer to these intra-subspace orbital ro- 
tations {EpE^, E£} as redundant, since the total wavefunc- 
tion is invariant under such orbital transformations, if accom- 
panied by the corresponding transformation of CI coefficients 



The redundant orbital rotations can be excluded in vary- 
ing our action functional in Eq. ( p"8| ), since their effects to 
6*i> are taken into account by the CI variations 5Cj. On 
the other hand, for {/i, v) belonging to different orbital sub- 
spaces (core-active, core-virtual, or active-virtual), the projec- 
tor IIcas can be dropped in Eq. (21 1, and we have a simpler 
expression, 



<*l 



H-R, E£ 



l*) = 0, 



(22) 



with E£ = {EL 1 , E%, E\,E\}, constituting the non-redundant 
orbital rotations. The general orbital-EOM of Eq. (21 1 is thus 
reduced to Eq. ( |22] >, which is to be solved only for the non- 
redundant orbital pairs. 



It is fascinating to see an analogy in Equation (22i with 



the time-independent MCSCF theory; It is formally identi- 
cal to the generalized Brillouin condition of the stationary 
wavefunction [73, 74|, if we replace H — R with H. Thus, 
the remaining derivations are parallel to the time-independent 
theory. We can explicitly write down the matrix elements of 



Eq. ( 22 1 to obtain 
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(24) 
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and P^ ee 



'ME"X\^>) are one- 



y pAl 



where D£ = 

and two-electron reduced density matrix (RDM) elements, re- 
spectively. The matrix F is called the generalized Fock ma- 
trix, whose Hermiticity, leading to vanishing right-hand side 
of Eq. ( |23] l, is the stationary condition with respect to the or- 
bital variations. g9]|5T]f53) 

The nonzero density matrix elements of the CASSCF wave- 



function are D) = 25),D l u ,P l k ] = 



4^ 



2SiSl,P% = 



2D*,P/^ — —D l u and P*™. Then the required generalized 
Fock matrix elements read 1 5 1 1 



Ft 



Ft 



2(/r+Gn, 



%d? + (r t y t 



(25) 



where the matrices /, G, and r t represent, respectively, oper- 
ators /, G, and ft given by 



f = h 



FC 



d.c. 



2Jj - Kj i . 



o=\Ji-\ki)D; 



r t \4> t ) = WZ\cf> v )P™, 



f.c. 

h FC {t) = h{t) + J2 (2^(0) - Kj(0) 



(27) 



(28) 



(29) 



(30) 



where summation j in Eqs. ( 27 1 and ( 30 1 are restricted within 



dynamical-core (d.c.) and frozen-core (f.c.) subspaces, re- 
spectively. The operators / and G are universal and Hermi- 
tian, while P t is defined with an active orbital <j> t to be ap- 
plied from the left, and non-Hermitian. We define Coulomb 
J, exchange K, and general W mean field operators as J p — 
JP , K p = kv , Jp = WP ,KP\<p r ) = WP \cj> q ) 
local [63 1 and given in the coordinate space as 



where W]! is 



WP(r)= df<j>;{f)V ee {r,f)<p q {f). 



(31) 



In Eq. pO) , time argument t is explicitly attached to empha- 
size that the time-dependence of the frozen-core dressed one- 
electron Hamiltonian h FC (£) comes entirely from the external 
laser field contribution in h(t). Now Eq. (23 1 for the time 
derivative matrix Rg = i(0 p |</>„) can be worked out for inter- 
subspace (non-redundant) elements: 



Rt = R l „ 



(ief.c), 



R? = B£ = f? + G* (ied.c), 



Rt = K* = ft + {r u )l (d- 



(32) 



(33) 



(34) 



TDt TDl 

K i - Rt 

with Dl = 251 
elements: 



(t G d.c), ^ ' 

-D^, and for intra-subspace (redundant) 



RH = (^\9(t)\4>„ 



(36) 



where 0(t) can be an arbitrary one-electron Hermitian oper- 
ator [38 59, 66 1, reflecting the invariance of the total wave- 
function against the redundant orbital transformations. 

One could, in principle, directly work with Eqs. d32|-(|36]) 
in the matrix formulation [59|, which determines the time de- 
pendence of occupied {0 p (£)}, as well as virtual {(f> a (t)} or- 
bitals. However, it is beneficial to introduce the orbital pro- 
(26) jector Q = J2 a \4>a){(t>a\ = 1 - J2 P I^pX^pI onto the virtual 



orbital space, to avoid (using the assumed completeness) ex- 
plicitly dealing with numerous virtual orbitals Il38ll67l . Thus 
we arrive at the final expression of EOMs for dynamical-core 
and active orbitals as follows: 



t)=Q {/>*) + ?»\K) (D- l ) u t ] + \<j> p )R p t , 



(37) 



(38) 



and R is determined by Eqs. (35 1 and (36i with a particular 
choice of 6(t). Solving these equations guarantees the opti- 
mal separation, in the TDVP sense, of frozen-core, dynamical- 
core, active, and virtual orbital subspaces, as illustrated in 
Fig.[T] This ensures the gauge-invariance of the TD-CASSCF 
method, since the orbital subspaces are stable against single 
excitations [Eq. (22 1] arising with the transformation, e.g., 



from the length gauge to the velocity gauge. 



2. CI equations of motion 



The general CI-EOM of Eq. ( 20 1 is specialized to the TD- 
CASSCF method as 



ricAs 



id j = J2 ( H u - s " eA - R ") °j> 



(39) 



where Rjj — (I\R\J), E A and Hfj are active orbital contri- 
butions to the total energy and determinant basis Hamiltonian 
matrix elements, respectively, 



e = miim = e c + e 



rA 



(40) 



D. Computational remarks 

The TD-CASSCF method includes as special cases both 
TDHF and MCTDHF methods, thus bridges the gap between 
the uncorrelated and fully correlated descriptions in a flexible 
way. A practical advantage of this generality is that a com- 
putational code written for the TD-CASSCF method can be 
used also for single-determinant TDHF and MCTDHF cal- 
culations, by setting {uq — n, n A = 0} and {tiq = 0, 
n A = n}, respectively. It can also execute open-shell TDHF 
calculation with fixed CI coefficients [66 1 . One indeed finds 
close similarity between TD-CASSCF EOMs Eqs. ( [37f|39| > 
and those of the MCTDHF method (See e.g. Ref. [65]). 
Naively, ingredients of the TD-CASSCF EOMs are the com- 
pilation of those for TDHF (core orbitals) and MCTDHF (ac- 
tive orbitals and CI coefficients) methods. This means that an 
existing code for the MCTDHF method can be easily general- 
ized to the TD-CASSCF method. 

The computationally most demanding procedures required 
to integrate the TD-CASSCF EOMs are grouped into two cat- 
egories: 

(A) Calculations of 2RDM elements P^, and the two- 



electron contributions of Eq. ( 39 1 



iCj 



1 



IIcAS 

2^ 

J 



9uw(^U)t 



'J f~l 

0.7 



(45) 



The amount of work in these procedures roughly scales 
as 0(N'l(n A - N A ) 2 N dct ) if N£ = JVf , (see Ref. E2 
for more details), where /Vdct is me number of determi- 
nants in ncAS which in turn scales factorially with the 
number of active electrons N A . 

(B) Calculations of the mean fields WjSjir), two-electron in- 



tegrals g™ w , and the 2RDM contributions in Eq. (38 1, 



(I\H\J)=6 U E C + Ht 



J' 



where 



£ c =e{(^ fc ):+^} 



pA ft r\u 1 1 tv puw 



H i j = fl( D u)t + o9uw(Pu)t™> 



(41) 



(42) 



(43) 



(44) 
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(46) 



whereCD/j)? = (I\El\J), and (P u )™ = {I\E%\J). In 
Eq. (39 1, we make a particular phase choice, i(^ r |\i r ) = 0, 

by extracting the dynamical phase cxp — ij dt'E(t') from 

the total wavefunction. This stabilizes the CI-EOM especially 
when we have a large active space. 



The computational cost of these steps depend explicitly 
on the number of grid points (or basis functions) Nb, 
as 0(n 2 N 2 ) for the mean fields and 0(n A iV{,) for the 
others. 

Important cost reductions are achieved for both procedures 
(A) and (B) by the TD-CASSCF method adopting core or- 
bitals, compared to the MCTDHF method with the same num- 
ber of occupied orbitals. The speed-up and resource savings 
for procedure (A) is substantial due to the decreased CI di- 
mension. This is especially the case if N A <C N, which is 
expected for an electronic structure with a few weakly-bound 
valence and large numbers of physically inactive core elec- 
trons. The cost reduction for procedure (B) is not as drastic 
as for (A), since the amount of arithmetics 0(n 2 N 2 ) of com- 
puting mean fields W?(r) is independent of the CAS structure 
(only related to n). The computations of two electron integrals 



and Eq. (46 » become much faster through restricting the or- 
bital indices within the active instead of all occupied orbitals. 



Relative importance of these bottlenecks largely depends on 
the problem at hand, and on the spatial representation of the 
orbitals and electron-electron interactions. This point will be 
discussed in Sec. lIIICl 



III. NUMERICAL RESULTS AND DISCUSSIONS 

In this section, we apply the TD-CASSCF method to the 
ionization dynamics of one-dimensional (ID) multielectron 
model molecules. The effective ID Hamiltonian for N elec- 
trons in the potential of M fixed nuclei interacting with an 
external laser electric field E(t) is taken as 



H 



N 

E 



i o 2 

'2dx? 



M 

E 



v^ 



x„ 



N 



1 



E(t)^\+J2^7 r 

J T^j y/iXi-Xj) 2 + d 



(47) 



where Xi (i — 1, • • • , N) is the position of the z-th electron, 
X = {X a } and Z = {Z a } (a = 1, • • • , M) are the positions 
and charges of nuclei, and c and d adjust the soft Coulomb op- 
erators of electron-nuclear and electron-electron interactions, 
respectively. The electron-laser interaction is included within 
the dipole approximation and in the length gauge. Note that 
the TD-CASSCF method is gauge-invariant as mentioned in 
Sec. ~" 



IIC 1 



We have performed some of the calculations de- 
scribed below also in the velocity gauge, and confirmed that 
the results are virtually identical to those in the length gauge. 
In this work, we make the simplest choice of 6(t) = 
in Eqs. (36 1, and therefore Ru = in Eq. (1391. The 



orbital-EOMs are discretized on an equidistant grid of spac- 
ing Ax — 0.4 (finer grid with Ax — 0.1) is used for drawing 
Figs. 2 4 1, within a simulation box \x\ < L with L = 600. An 



absorbing boundary is implemented by the mask function of 
cos 1 / 4 shape at 15% side edges of the box. The ground-state 
electronic structure is obtained by the imaginary time prop- 
agation with the fourth-order Runge-Kutta (RK4) algorithm 
with Schmidt orthonormalization of orbitals after each prop- 
agation |59|. The real-time propagations use variable step- 
size embedded fourth- and fifth-order Runge-Kutta (VRK5) 
method. The kinetic energy operator — | J^ is evaluated by 

the eighth-order finite difference, and spatial integrations are 
replaced by grid summations using the trapezoidal rule. Fur- 
ther details of the computations are given separately below. 



A. lD-LiH and LiH dimer models: Ground-state 

We consider ID lithium hydride (LiH) and LiH dimer mod- 
els. The reason for choosing these models is that they rep- 
resent the simplest examples of such electronic structures 
with (i) deeply bound orbitals, and (ii) several weakly bound 
orbitals, as shown below. These characteristics should be 
the key in the three-dimensional (3D) multielectron dynam- 
ics, where the existence of energetically closely-lying valence 
electrons is quite common, which requires to take both multi- 
channel and multielectron effects into account. As discussed 
previously [62|, cares have to be made for the physical sound- 
ness of ID models. Nevertheless, we expect that the features 
(i) and (ii) are transferable, and ID applications can elucidate 
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FIG. 2. Several cuts of the adiabatic energy surface E(r, R) of 
the lD-(LiH) 2 model obtained by the MCTDHF(8) method. The 
total energy of Eq. {481 is plotted against intermolecular LiH- • LiH 
distance R for several bond lengths r of LiH, constrained to be the 
same for two LiH molecules. 



advantages and limitations of theoretical methods, before ap- 
plied to real 3D systems. 

For LiH, we set molecular parameters as Z = {3, 1} and 
X = {-1.15,1.15}. For (LiH) 2 , Z = {3,1,3,1} and 
X = {-4.05, -1.75, +1.75, +4.05}. The soft Coulomb pa- 
rameters c = 0.5 and d — 1 are used, since an often-made 
choice of c = d = 1 1751 l76l was found to overemphasize the 
electron-electron repulsion. The above molecular parameters 
correspond to the equilibrium bond length r = 2.3 of Li-H 
and intermolecular distance R = 3.5 of LiH-LiH, as shown 
in Fig. [2] which plots several cuts of the adiabatic energy sur- 
face of (LiH>2, 



M 



E(r,R) = {$\H\*) + J2 



a>b 



Z a Zt 

X„ — Xh 



(48) 



The energy surface with parameters c = d = 1 predicted no 
stable LiH dimer in this nuclear configuration, relative to the 
separated LiH molecules. 

To make a sensible comparison among methods with dif- 



ferent active spaces, we consider the following wavefunctions 
for LiH: 



HF: core ■ 

HF: valence ■ 

HF: total - 

MCTDHF(8): total ■ 



*hf : fyl, 


(49a) 


*CASSCF(2,n) : 4>\ (0203 — 4>n+l) , 


(49b) 


^MCTDHFfn+l) = (010203— 0n+l) , 


(49c) 


and for (LiH) 2 : 




* HF : 4>\4>l4>l4>l, 


(50a) 


*CASSCF(2,n-l) = 0?0203 (0405-0n+2) , 


(50b) 


^CASSCF(4,n) ■ 0?02 (0304-0n+2) , 


(50c) 


*MCTDHF(n+2) ■ (010203 — 071+2) , 


(50d) 



following the notations of Eqs. ^ and Q. The CASSCF and 
MCTDHF wavefunctions are designed to consist of the same 
number of occupied orbitals with increasing active orbitals. 

Figures [5]and|4]show the shapes of the ground-state occu- 
pied HF orbitals and the one-electron probability density for 
LiH and (LiH)2, respectively. As seen in Fig. |3lb) the node- 
less first deepest HF orbital of LiH localizes at Li "atom", 
while the second orbital is responsible for the formation of 
"chemical bond", made from constructive superposition of the 
ground-state wavefunction of H and the second atomic orbital 
of Li, the node of the latter shifted to the bonding region. 
Figure [3|a) shows that the total electron density is well re- 
produced by HF method compared to the MCTDHF(8) den- 
sity. In Fig. Kb), one sees that HF orbitals of (LiH)2 can be 
clearly separated to the deeply-bound core (orbitals 1 and 2) 
and weakly-bound valence (orbitals 3 and 4) orbitals, the for- 
mer keeping the atomic-orbital characters of Li, while the lat- 
ter two orbitals delocalizing across the dimer. Again, as seen 
in Fig.Q a), the total density is well reproduced by HF. Finally, 
it is observed that the tails of the total electron density are de- 
termined by the valence electrons both in LiH and (LiH)2. 

Table II] summarizes the ground-state calculations. As seen 
in the table, there are significant gaps in total energies between 
methods with different numbers of active electrons. However 
the MCTDHF values for the other properties are reproduced 
rather well, by CASSCF(2, n) and CASSCF(4, n) methods for 
LiH and (LiH) 2 , respectively. For instance, the difference in 
JStot of CASSCF(4, 8) and MCTDHF(IO) for (LiH) 2 is ap- 
proximately 12 mHartree, but those in IP and A 2 are 0.05 eV 
(2 mHartree) and 0.003 eV (0.2 mHartree), respectively. This 
tells that the correlations responsible for these properties are 
those among the valence electrons. 

The CASSCF(2, n) active spaces of Eq. pOb} , with only 
two of four nearly degenerate electrons being correlated, are 
not physically sensible ones for (LiH) 2 . Accordingly, the re- 
sulting dipole moment values are not much improved from 
the HF value. More seriously, the proper dissociation limit 
of such wavefunction to the equivalent LiH molecules can- 
not be defined well, i.e., the formation energy A 2 cannot 
be obtained. This problem is due to the lack of the "size- 
extensivity", ll48l l49l l53ll which fails to guarantee the equal 
quality of the approximation for different electronic config- 
urations. The size-inextensive treatment covers less and less 
electron correlation as systems grow larger. 
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FIG. 3. The ground-state electron density and occupied HF orbitals 
of the lD-LiH model, (a) The total electron density by HF (black 
dashed) and MCTDHF(8) (black thick solid) methods are compared. 
The two lines closely overlap each other. Also shown are the core 
(orbital 1) and valence (orbital 2) contributions to the total HF den- 
sity, (b) The occupied HF orbitals are drawn in an arbitrary scale 
vertically shifted by orbital energies -1.82 (orbital 1) and -0.67 (or- 
bitals 2). The solid curve and dashed vertical lines show the nuclear 
potential and positions of nuclei, respectively. 
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FIG. 4. Same as Fig. [3] for the lD-(LiH)2 model, (a) The core 
and valence densities correspond to the contributions from orbitals 1 
and 2, and orbitals 3 and 4, respectively, (b) The orbital energies are 
-1.85, -1.77, -0.73, and -0.60, for orbitals 1 to 4, respectively. 



B. ID LiH model: Ionization dynamics 



Now we apply the TD-CASSCF method to the laser-driven 
electron dynamics of the lD-LiH model. We use the three- 



cycle laser electric field of the following form; 

t 



E(t) = E sin(wt) sin 2 I tt- ] , < t < r (51) 

with u) = 0.06075 (wavelength 750 nm), r = 6ir/ui, and 
three different amplitudes E = 0.0534, 0.107, and 0.151 cor- 
responding to peak intensities Iq = l.Ox 10 14 , 4.0x 10 14 , and 
8. Ox 10 14 W/cm 2 , respectively. The Keldysh parameters are 
1.30, 0.65, and 0.46, respectively, for the three intensities. In 
view of the ground-state electronic structure of Fig. [3] and the 
above laser profile, one reasonably expects that the dominant 
physical process involved is the tunneling ionization from 
the highest occupied orbital in the static HF picture. Hence, 
we can speculate that the two-active-electron description TD- 
CASSCF(2, Ha) is necessary and sufficient for the accurate 
description of the dynamics, as will be confirmed below. 
Figure E\ shows the time evolution of the dipole moment. 



TABLE I. Total energy E to t in atomic unit (a.u.), dipole moment (x) 
in a.u., first ionization potential IP in eV, LiH formation energy Ai in 
eV, and (LiH)2 formation energy A2 in eV. Results of HF, CASSCF, 
and MCTDHF methods with varying active spaces are compared. 
The ionization potential (IP) is computed as the difference of the to- 
tal energies of the neutral and cationic ground-states. The Ai is the 
difference of the energy of LiH and the sum of energies of Li and H, 
and A2 is the difference of the energy of (LiH)2 and twice that of 
LiH. For computing Ai and A 2 with CASSCF and MCTDHF meth- 
ods, the active spaces for the fragments are the proper dissociation 
limit of the parent description of the complex. For HF calculations, 
the open-shell restricted HF method is used for doublet species. The 
A2 of CASSCF(2, ha) methods are not available since their proper 
dissociation limits are not well defined. 
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First we observe the large difference in the results of TDHF 
and other methods. For the lowest intensity of l.OxlO 14 
W/cm 2 , the difference remains quantitative, largely due to 
the difference of the ground-state permanent dipole moment. 
For higher intensities, TDHF clearly underestimates the laser- 
driven large-amplitude electron motions. This is due to 
the fundamental inadequacy of the closed-shell description, 



Eq. (49a 1, of the tunneling ionization process, which involves 
spatially different motions of the ionizing and non-ionizing 
electrons. The TD-CASSCF(2, 2) brings substantial improve- 
ment over the TDHF, giving results with much better agree- 
ment with the MCTDHF ones. The convergent description in 
the TD-CASSCF(2, n A ) series is obtained at n A = 4. The 
TD-CASSCF(2, tia) with n\ > 4 closely reproduce the re- 
sults of MCTDHF method. 

Figure |6]plots the n-electron ionization probability P n , de- 
fined for convenience as a probability to find n electrons lo- 
cated outside a given distance R lon — 20 (see Appendix [A|, of 
LiH as a function of time for the peak intensities (a) 4x 10 14 
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FIG. 5. Dipole moment of the lD-LiH model as a function of 
time, with peak intensities (a) 1 x 10 14 , (b) 4x 10 14 , and (c) 8x 10 14 
W/cm 2 . 
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and (b) 8x 10 14 W/cm 2 . No appreciable ionization is found 
with the lowest intensity. The probability of finding more than 
two ionized electrons is negligibly small for all intensities. As 
seen in Fig.|6| TD-CASSCF(2, 4) gives virtually the same re- 
sults as MCTDHF(5). The TDHF method, on the other hand, 
underestimates single ionization Pi and, at the higher inten- 
sity, unphysically overestimates double ionization P 2 . This 
is the consequence of forcing two valence electrons to travel 
with a single spatial orbital. 

These results demonstrate that the TD-CASSCF(2, 2) con- 
stitutes the simplest method to describe the present dynamics 
in a physically correct way. Its total wavefunction can be writ- 
ten as 



W = A [0i0i {Ci0 2 02 + C 2 0303}] 



(52) 



in the natural orbital representation 17711781 . where <j>i (<fii) is 
an orbital occupied by up (down) spin electrons. The two- 



configuration CI part of Eq. ( 52 1 can be transformed back to 
the non-orthogonal expression [Eq. ([Ta|], 



* oc A 



' i " i <| ("02 ^3 + ^3 ^2) —p- («/3 - Pa) 



(53) 
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FIG. 6. Ionization probabilities P n of the lD-LiH model as a 
function of time, with peak intensities (a) 4x 10 14 and (b) 8x 10 14 
W/cm 2 . TDHF (black dotted), TD-CASSCF(2, 4) (colored), and 
MCTDHF(5) (block solid) results are compared. 



giving a clearer picture of different spatial motions of the two 
valence electrons, with [48] 



K^ihMI 



|Ci|-|C 2 | 
ICil + ICal 



(54) 



The flexibility inherent in Eqs. (52 1 or (53 1 enables a seam- 
less transition from the closed-shell dominant ground-state 
(|Ci| 3> | C 2 1 •<=>■ (ipi\ip2) ~ 1) to the single ionization limit 
(|Ci| ~ |C 2 | <=>■ ("0i|02) ~ 0). The ionization dynamics, 
therefore, is characterized by the strong or static correlation 
B8ll49ll52ll53l in the sense that it involves drastic changes of 
the configuration weights (the magnitudes of CI coefficients) 
with more than one determinants contributing significantly. 
The failure of single-determinant TDHF to describe the ion- 
ization process is attributed to the lack of this type of correla- 
tion. 

For quantitatively accurate description of the dynamics, the 
above minimum CI wavefunction has to be improved by incor- 
porating more-than-two active orbitals, as seen in the conver- 
gence of the dipole moments in Fig.|5]with respect to the num- 
ber of active orbitals. The agreement of TD-CASSCF(2, ua) 
and MCTDHF results indicates that the core electron corre- 
lation is not relevant, at the first approximation, for the ion- 
ization dynamics induced by the present laser field. The TD- 
CASSCF allows the compact representation of such physical 
situations. 



C. lD-LiH dimer model: Ionization dynamics 

In this section, we proceed to the multielectron dynamics 
of lD-(LiH) 2 model. We assess TDHF, TD-CASSCF(2, 7), 
TD-CASSCF(4, 8), and MCTDHF(IO) methods. These active 
spaces are shown in Eq. ( |50| ) with n = 8. The latter two are 
twice the size of those in TD-CASSCF(2, 4) and MCTDHF(5) 
for LiH, respectively, which have been confirmed to provide 



the convergent description in Sec. IIIB 



Figure I7J shows the temporal evolution of the dipole mo- 
ment simulated with various methods. One clearly sees that 
TDHF and TD-CASSCF(2, 7) results show large deviations 
from MCTDHF(IO) ones, while TD-CASSCF(4, 8) repro- 
duces the results of MCTDHF(IO) fairly well. This indicates 
that all the four valence electrons sketched in Fig. |4] actively 
participate in the field-induced ionization dynamics (this does 
not necessarily mean that the four electrons are ionized), while 
tightly bound core electrons remain non-ionized. For the ion- 
izing electrons, the closed-shell description is inadequate as 
discussed in Sec. IIIIBI 

In Figs. [8] and [9] we compare the temporal evolution of the 
ionization probability P n with i?; on = 20 of (LiH) 2 com- 
puted by approximate methods and MCTDHF(IO). As can be 
seen from Fig. [8] both TDHF and TD-CASSCF(2, 7) methods 
tend to underestimate single ionization for all the examined 
intensities. The probability of finding more than two ionized 
electrons are found to be erroneous in an inconsistent way, 
thus not shown. In a striking contrast, TD-CASSCF(4, 8) 
well reproduces the ionization probability P n obtained with 
MCTDHF(IO) [Fig|9](a)-(c)]. Slight deviation is seen only 
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FIG. 7. Dipole moment of the lD-(LiH)2 model as a function of 
time, with peak intensities (a) 1 x 10 14 , (b) 4x 10 14 , and (c) 8x 10 14 
W/cm 2 . Results of TDHF, TD-CASSCF(2, 7), TD-CASSCF(4, 8), 
MCTDHF(IO) methods are compared. 



at the later stage of the pulse for the higher intensities. The 
inclusion of more active orbitals would further improve the 
agreement. 



So far, all the core orbitals have been treated as dynami- 
cal core. In Fig.|9](d)-(f), the ionization probability computed 
with TD-CASSCF(4, 8) with all the core orbitals treated as 
frozen, denoted TD-CASSCF(4, 8)-FC, are shown. It repro- 
duces the results of MCTDHF(IO) almost as nicely as TD- 
CASSCF(4, 8) with dynamical-core orbitals, which indicates 
that the core polarization plays minor roles in the present dy- 
namics. 



It is worth noting that even at the lowest intensity 
l.OxlO 14 W/cm 2 dominated by single ionization, the TD- 
CASSCF(2, 7) fails to give an accurate value of P\ but un- 
derestimates it roughly by half [Fig. [81(d)]. This implies the 
importance of the multichannel ionization, which can be de- 
scribed correctly only when all the relevant orbitals are in- 
cluded in the active space. On the other hand, at higher inten- 
sities, the total wavefunction consists of the widespread super- 
position of the ground-, excited- and continuum-states. For 
a balanced description, each of these components has to be 
treated with an equal quality, which requires a size-extensive 
theory. The MCTDHF, as the exact theory within a given 
number of time-dependent bases, fulfills the size-extensivity 
condition. The TD-CASSCF with a proper active space pre- 
serves this important property of the MCTDHF. It is demon- 
strated by the accurate multiple ionization probabilities ob- 
tained by the TD-CASSCF(4, 8) method, up to P 4 for the 
highest intensity in Fig. [9] (c). The importance of selecting 
an appropriate active space is illustrated by the fact that the 
TD-CASSCF(4, n A ) is required for (LiH) 2 , while the TD- 
CASSCF(2, n A ) is adequate for LiH. 



D. Analyses of computational cost 



Table |TT] summarizes computational times for simulations 
of the lD-(LiH) 2 model with a peak intensity 4x 10 14 W/cm 2 . 
To highlight the different computational bottlenecks discussed 
in Sec.|lfc] several box sizes (N b = 1000, 2000, and 3000) 
are considered. The CPU times in table |llj-(a), (b), and (c) are 
recorded on a single Xeon processor of clock frequency 3.33 
GHz for propagating 1000 time-steps during 2T < t < 2. IT 
with the fixed step-size RK4 algorithm, where T — 2tt/u>. 
Entry (d) compares wall clock times spent for completing the 
simulation up to four optical cycles (0 < t < AT), with the 



r 



VRK5 algorithm, measured for multi-threaded computations 
using 12 processors. 

First, as seen in table |ll[(a), CPU times for procedure 
(A) grows rapidly with increasing CI dimension, reproduc- 
ing the theoretical linear dependence on A^jet- These timings 
marginally depend on Nb. Next, CPU times for procedure (B) 
in table|n[(b) scale as O(N^) with d = 1.95, 1.67, 1.66, and 
1.47 for TDHF, TD-CASSCF(2, 7), TD-CASSCF(4, 8), and 
MCTDHF(IO) methods, respectively. This is the consequence 
of competing 0(n 2 N£) and 0(n\Ni>) contributions as dis- 
cussed in Sec. |IID| with growing importance of the latter for 
larger active spaces. The TD-CASSCF(4, 8)-FC demands less 
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FIG. 8. Ionization probabilities Pq and Pi of the lD-(LiH)2 model as a function of time, with peak intensities lx 10 14 (left), 4x 10 14 (center), 
and 8x 10 14 W/cm 2 (right). Results of TDHF (top) and TD-CASSCF(2, 7) (bottom) are compared with those of MCTDHF (black solid and 
dotted). 
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FIG. 9. Ionization probabilities P n with n < 4 of the lD-(LiH)2 model as a function of time, with peak intensities 1 x 10 14 (left), 4x 10 14 
(center), and 8xl0 14 W/cm 2 (right). Results of TD-CASSCF(4, 8) (top) and TD-CASSCF(4, 8)-FC (bottom) are compared with those of 
MCTDHF (black solid). 



CPU times than the TD-CASSCF(4, 8), due to the strict lo- cality of frozen-core orbitals, limiting the range of exchange 
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operators Ki\4>t) around the core region. 

Net CPU times are listed in table |H}(c). In TDHF and TD- 
CASSCF calculations with core subspaces, the grid-intensive 
procedure (B) is definitely rate-limiting. In contrast, MCT- 
DHF calculations involve severe bottlenecks both in proce- 
dures (A) and (B). The iVdct -dependent works dominate 85%, 
67%, and 55% of the net CPU times, with N b = 1000, 2000, 
and 3000, respectively. The cost reduction achieved by the 
TD-CASSCF method largely depends on the relative impor- 
tance of procedures (A) and (B). Ratios of net CPU times for 
TD-CASSCF(4, 8) and MCTDHF(IO) calculations are 0.12, 
0.29, and 0.45 with N b = 1000, 2000, 3000, respectively. Sim- 
ilar trends are observed for wall clock times with VRK5 algo- 
rithm, as seen in table Inkd). The stability of EOMs is found 
to be similar for the tested methods, requiring about 70000 
evaluations of EOMs. The cost gain by the TD-CASSCF 
method relative to the MCTDHF method will be more dras- 
tic if Na <C N. However, an efficient implementation of the 
mean field potential [Eq. pT} ] is essential to achieve further 
speed-up for large N b , especially in three-dimensional appli- 
cations. 



TABLE II. Computational times for simulations of the lD-(LiH)2 
model. First row: Numbers of determinant TVdct within the symme- 
try of zero spin projection. Entries (a), (b), and (c): Central pro- 
cessor unit (CPU) times in minutes for (A), (B), and overall proce- 
dures denned in Sec. |IID[ respectively. Entry (d): Wall clock times 
in minutes spent to complete the propagation of four optical cycles. 
Different simulation boxes are employed: L — 200 (Nt, = 1000), 
L = 400 (N b = 2000), and L = 600 (N b = 3000). See text for 
more details. 
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N b = 3000 0.0 


0.3 


1.4 


1.4 


215.3 


(b) RK4 / 1000 steps, 


CPU-B 








N b = 1000 4.3 


26.8 


28.4 


22.0 


35.7 


N b = 2000 16.3 


87.6 


90.9 


73.8 


103.5 


N b = 3000 36.8 


168.0 


175.4 


132.9 


177.2 


(c) RK4 / 1000 steps, 


CPU net 








N b = 1000 4.4 


27.0 


29.8 


23.4 


248.9 


N b = 2000 16.4 


88.1 


92.5 


75.4 


319.4 


N b = 3000 36.9 


168.7 


177.3 


134.7 


393.5 


(d) VRK5 / 4 cycles, 


Wall 








N b = 1000 8.2 


43.7 


51.3 


39.1 


451.5 


N b = 2000 31.0 


174.8 


192.6 


141.4 


628.8 


N b = 3000 65.2 


378.5 


394.0 


282.9 


823.1 



IV. CONCLUSIONS 

We have developed a new ab initio time-dependent many- 
electron method called TD-CASSCF. It applies the concept 
of CASSCF, which has been developed for the electronic 
structure calculation in quantum chemistry, to the multielec- 



tron dynamics in intense laser fields, introducing frozen-core, 
dynamical-core, and active orbital subspaces. The classifica- 
tion into the subspaces can be done flexibly conforming to 
simulated physical situations and desired accuracy, and both 
TDHF and MCTDHF methods are included as special cases. 
This feature enables compact yet accurate representation of 
ionization dynamics in many-electron systems, and bridge the 
huge gap between TDHF and MCTDHF methods. 



We have applied the TD-CASSCF method to the ioniza- 
tion dynamics of lD-LiH and lD-(LiH)2, to assess its capa- 
bility to describe multichannel and multielectron ionization. 
It has been confirmed that the present method closely repro- 
duces rigorous MCTDHF results if the active orbital space is 
properly chosen to include appreciably ionizing electrons. We 
have also confirmed that the TD-CASSCF provides substan- 
tial computational cost reduction in the Cl-length dependent 
procedures, which scale by far the steepest with the system 
size in the MCTDHF method. Therefore, the TD-CASSCF 
method is most advantageous for problems in which only a 
few weakly-bound electrons out of a large number of total 
electrons ionize. 



While it is sometimes stated that the MCTDHF method is 
a time-dependent version of the CASSCF method ll60l I6D . 
this statement is even more suitable for the TD-CASSCF 
method introduced in the present study. With reduced compu- 
tational cost, the TD-CASSCF method with a properly cho- 
sen active space preserves most of the theoretically impor- 
tant properties of the MCTDHF: (i) flexibility to account for 
the strong-correlation involved in the ionization dynamics, (ii) 
size-extensivity, essential for a balanced description of differ- 
ent electronic configurations, (iii) gauge-invariance by virtue 
of the time-dependent variational optimization of orbitals, and 
(iv) invariance against orbital transformation within an orbital 
subspace, allowing e.g., the natural orbital analyses of the 
time-dependent wavefunction ||79l . 



It should be noted that the computational cost of the TD- 
CASSCF method still scales factorially with the number of 
active (not total) electrons, thus its applications are limited to, 
say, 16 half-filled active orbitals in view of the present state 
of the art in quantum chemistry. An example requiring such a 
large active space is the ionization from densely lying multi- 
ple valence orbitals in weakly-interacting molecular clusters. 
To approach to such a problem, more restricted (instead of 
complete) constructions of the active space will be necessary. 
Moreover, a breakthrough is needed to represent one-particle 
wavefunctions in the general molecular potential without par- 
ticular symmetries. In spite of these challenges, we foresee 
that the TD-CASSCF method will find fruitful applications in 
multielectron dynamics of, e.g., rare gas atoms heavier than 
helium, or molecules composed of atoms in the first few rows 
of the periodic table, exposed to visible-to-mid-infrared high- 
intensity pulses, which are inaccessible with the all-electron- 
active MCTDHF method. 
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etc, and Sfj is an TV x N matrix with its {ij} element being 
the inner-region overlap integral, 



(S< 



u)v 



dx<l>* p (i,I){x)<i>q{j,J){%) = (</>p|0g)<,(A5) 



Appendix A: Calculation of ionization probabilities 

To conveniently evaluate the multiple ionization yield in 
many electron systems, we introduce a domain-based ioniza- 
tion probability P n , defined as a probability to find n electrons 
in the outer region \r\ > i?i on and the remaining N — n elec- 
trons in the inner region \r\ < R[ on , with a given distance R 
from the origin, 



P 'n = ( J / dxi • • / dx n / dx n+1 
^*(xi,--,x N ) *^(Xi,-,X N ), 



dx 



N 



(Al) 



where J and J symbolize integrations over a spatial-spin 
variable x — {r, £ } with the spatial part restricted to the do- 
mains \r\ < R[ on , and \r\ > i?; on , respectively. 

It is convenient to introduce an auxiliary quantity T n ob- 
tained by replacing the outer-region integrals in Eq. \A\\ with 
the full-region ones (J — > J). It relates to P n as 



Pn 



E 

fc=0 



N ■ 



- n 
k 



{-l) k T n - k . 



(A2) 



By adopting the CI expansion of Eq. ([8]), and making use of 
the orthonormality of spin-orbitals in the full-space integra- 
tion, we have 



T n = -£ctCjD<8, 



(A3) 



/./ 



where 



D 



(o) 

ij 



N 

E 

N 

D IJ = /_^ € ij u " f ■'", I 
ij 

N N 
-)(2) _ X^S^JJJJ 



det (Sfj) , 

det (S7 



il) 



D iJ = E E ^Yi det (Sfj[ij :kl}), (A4) 



i>j k>l 



with 



MiJ) 



being the i-th (in a predefined order) spin-orbital 



in the determinant I. Sfj\ij ■ ■ : kl ■ ■] is the submatrix of S' < 



obtained after removing rows i, j, 
the latter, and 



and columns k, I, ■ 



u 
from 



JJ 



°qU,J) 



(-1) 4+J . 



(A6) 



The matrix Sfj and its submatrices are block-diagonal 
due to the spin-orthonormality, so that, e.g., det (Sfj) = 



j° 



I det 



det (S< 
determinant /. 



K^iPjn, 



e.g.. 
where I a is the a-spin part of the 



The procedure given above remains a manageable task in 
the present applications, up to eight (all) electron ionization 
probabilities in the lD-(LiH)2 model. While this scheme be- 
comes impractical for systems with more electrons, it may still 
be useful for problems where only a few electrons are ejected 
appreciably, since the dimension of Eqs. \K<\\ can be reduced 
to the number of the ionizing electrons. 

This approach allows the evaluation of multiple ioniza- 
tion yields by using the information of the inner region or- 
bitals (<f> p \(f>q)< and the formal orthonormality relation 6? = 
(4>p\4>q)< + {4 l p\ < t > q)>- ^ works with a reasonable size of the 
simulation box L, provided that i?i on <C L and a good ab- 
sorber is implemented to prevent the reflection of the wave- 
function. In fact, we performed calculations for the lD-(LiH) 2 



model in Sec. Ill using smaller boxes with L = 200 and 400 
a.u., where a sizable portion of the norm is lost at the bound- 
ary, and confirmed that the obtained ionization yields are vir- 
tually the same with those of Fig. [8] and [9] Such small- 
scale calculations could serve as preliminary validations for 
the choice of the active space before stepping into large-scale 
computations. 
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